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A stochastic system where bistability is caused by noise has been recently investigated by Bian- 
calani et al. (PRL 112:038101, 2014). They have computed the mean switching time for such a 
system using a continuous Fokker-Planck equation derived from the Taylor expansion of the Master 
equation to estimate the parameter of such a system from experiment. In this article, we provide 
the exact solution for the full discrete system without resorting to continuous approximation and 
obtain the expression for the mean switching time. We further extend this investigation by solving 
exactly the Master equation and obtaining the expression of other quantities of interests such as the 
dynamics of the moments and the equilibrium time. 


I. INTRODUCTION. 

In some stochastic systems, noise can have counter in¬ 
tuitive effects and the behavior of the system be markedly 
different from its deterministic, mean field approxima¬ 
tions. In some oscillatory gene networks, the regular os¬ 
cillations are caused by noise and cease in their absence 
[1]. In population genetics, the noise term can explain the 
emergence of less fit “altruistic” individuals Q. In ecol¬ 
ogy, the spatial aggregation of individuals can be caused 
by noise |3j, ij ; a similar explanation lies behind neutron 
clustering in nuclear reactors [5]. 

The general theory of noise induced transition in non¬ 
equilibrium systems has been extensively investigated by 
Horsthemke and Lefeve [bj. In the context of chemical 
equations and specifically genetic regulatory networks, 
there has been an intense investigation of systems where 
bistability is caused by noise and is absent from the deter¬ 
ministic formulation of kinetic rate equations. Samoilov 
et al. @] have considered the enzymatic futile cycle re¬ 
action and have shown that addition of noise can cause 
bistability and dynamic switching in the concentration of 
the substrate. Artyomov et al. [8J have considered a sim¬ 
ple model of T cells response and have shown again that 
in the presence of noise, the steady state distribution can 
become bi-modal. Qian et al. §}] and Thomas et al. 0, 
using different approaches, have derived a general frame¬ 
work to elicit the role of fluctuation time scales separa¬ 
tion in the appearance of noise induced bistability. In an 
elegant experiment, To and Maheshri 0 have investi¬ 
gated a synthetic transcriptional feedback loop and have 
demonstrated the bimodality of the response without co¬ 
operative binding of the transcription factor, a usual hy¬ 
pothesis to explain bistability of genetic switches. 

Recently, Biancalani et al. [12j investigated another 
stochastic system where bistability is caused by noise: in 
this system, individuals (or molecules) can be in one of 
the two configurations A and B and can switch from one 
to the other according to the following transition rates: 

W~(n) = W(n —> n — 1) = [r(N — n) + e) n (1) 

W + {n) = W(n —»• n + 1) = (rn + e) (N — n) (2) 

where n is the number of individuals in configuration A 


and N is the total number of individuals. In the follow¬ 
ing, n is used to characterize the state of the stochastic 
system at a given time. The rate r characterizes the two 
body interactions 

Xi + Xj A 2X t i = A, B- j = B,A 

while the rate e characterizes spontaneous switching of 
an individual from one configuration to the other: 

Xi —> Xj 

Without loss of generality, we will set r = 1 in the fol¬ 
lowing. This is achieved by scaling both time and e by 
the factor r. 

Such a system can model for example a colony of forag¬ 
ing ants collecting food from two sources. In population 
genetics, this is the Moran model for two competing al¬ 
leles A and B with bidirectional mutations [H|. Such 
systems were also proposed in the context of autocat- 
alytic chemical reactions with small number of molecules 
EUm , or the dynamic Ising model 0 for a set of fully 
connected spins. The general properties of this stochas¬ 
tic system, and its application to population genetics in 
fluctuating environment were discussed by Horsthemke 
and Lefeve [6j. 

The behavior of this system is markedly different from 
its mean field, deterministic approximation. Indeed, the 
equation for (n) , the mean number of individuals in one 
state, is: 

^± = (W + (n)-W-(n)) = e(N-2(n)) (3) 

and has a stable stationary solution (n) = N/2. However, 
for small values of e, i.e. e -C 1 /N, the system is observed 
most of time in one the two boundary states n = 0 or 
n = N, and seldom in states close to n = N/2. The 
bistability of the system is caused solely by the noise and 
cannot be captured by the mean field equation ([3]) . 

The reason behind the bistability is the following: 
in the absence of spontaneous switching (e = 0), the 
states n = 0 (all individuals in configuration B ) and 
n = N (all individuals in configuration A) are absorb¬ 
ing: W + (0) = W~(N) = 0. Eventually, the sys¬ 
tem will end up in one of these two states and remain 
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there. When e > 0, these states cease to be absorbing. 
However, the mean residence time r in these states is 
(W + (a) + W~{a)) 1 = 1/eN (where a = 0, AT) while 
the residence time in other states is 0(1). Therefore, in 
the regime eN <C 1, the system is observed mostly in the 
boundary states. 

In their article, Biancalani et al. computed T(0), the 
mean switching time (the mean first passage time) from 
state n = 0 to state n = N, and show that the observa¬ 
tion of this quantity can lead to the measure of the pa¬ 
rameter e of this stochastic system. For this computation, 
they expanded the Master equation of the stochastic sys¬ 
tem in powers of 1/IV and neglected terms of 0(1 /TV 3 ) to 
obtain the forward and backward Fokker-Plank equation, 
from which the mean switching time can be obtained ( 
[nj , equation (4) and Supplementary Materials, equa¬ 
tions (4) and (11) ). This approximation is fragile, spe¬ 
cially for small N where the noise is strong. In particular, 
to compute T(0), they have used two different approxi¬ 
mations, one of which is valid for 0.2 < Ne and the other 
for Ne —> 0, and there is no clear criterion for their over¬ 
lap. In this article, we compute the exact expression for 
T( 0) without any approximation, which is valid for all 
values of e. We further extend this investigation by giv¬ 
ing the exact solution of the discrete Master equation 
through the use of the probability generating function 
associated to the probabilities. Other quantities that we 
compute, such as the dynamics of the moments or the 
dynamics of the boundary states probabilities, provide 
other useful tools to measure and investigate this system. 

This article is organized as follow: in the next section, 
we give the exact expression for the mean first passage 
time T(n). The following section is devoted to the solu¬ 
tion of the Master equation. The final section is devoted 
to discussion and conclusion. 


II. SWITCHING TIME. 


equation for t{x ) which can be solved in terms of the 
hypergeometric function, as was done by Biancalani et 
aljl2l| fsee IV Bl) . The continuous limit is however fragile 
when e —>• 0, and the first solution obtained by Biancalani 
et al. does not converge to the right value in this limit. 
This is due to the absorbing boundary condition f'(0) = 0 
used in the continuous approximation, which fails in the 
limit e —> 0 as it can be observed directly from equation 
(ED (see also [12] Supplement. Materials). In order to 
resolve this problem, they have resorted to a limit process 
for the case e —> 0 by approximating (53 Supplement. 
Materials, eq.(28) ) 


2-Fi(^;|; lH ! 2e ) «2i i i(i,u;|;l) 


where u = Ne or 1 — Ne, i.e. setting e = 0 in the fourth 
argument of the hypergeometric function, but not in the 
second. This ad hoc approximation gives the correct so¬ 
lution for e —>• 0; no criterion however can be obtained 
for the overlap between the two solutions (figure [2D ■ 
These complications are due to the continuous approx¬ 
imation and can be avoided if the solution is computed 
directly for the discrete equations m■ The discrete so¬ 
lution is computationally much simpler, is valid for the 
whole range of e and N and does not involve any approxi¬ 
mation; specifically, the boundary conditions are set nat¬ 
urally and don’t need to be adjusted as a function of e. 
The solution is obtained by setting yj. = T(k) — T(k — 1), 
which transforms equations m into a simple one-term 
recurrence equation. The exact solution is then 


Vk+i = 


E 


2=0 


(N — k + e)(k-i) (i + l)(fc-j) 
(N — fc)(fc_i+ 1 ) (« + e)(k-i+i) 


0 <k<N 


where (a)( m ) = a(a + l)...(a + m — 1) = F(a + m)/T(a) 
is the Pochhammer symbol. 

As T(N) = 0, the first passage times T{in) are easily 
recovered from the yi~ : 


Preparing the system at time t = 0 in the initial state 
n = m, the system evolves and will reach the state n = N 
for the first time at some time T(m). The mean first 
passage times T[m) are obtained from the backward Kol¬ 
mogorov equation and form the linear system [18] 

W+(0) (T(1)-T(0)) = -1 (4) 

W + (m)(f(m + l)-T(m)) + 

W~(m) (f(m-l)-T(m)) = -1 (5) 

where 0 < m < N. Note that as W~( 0) = 0, we don’t 
need to write a separate equation ED f° r the boundary 
term T(0) ; the above notation however is clearer and 
highlights the boundary condition. Note also that by 
definition, T(N) = 0, so the above square system of linear 
equations is well posed. 

Using the continuous approximation n —>• x = n/N , 
T{m) —> t(x), and developing equation ED to the second 
order in (1/iV), one obtains the second order differential 


N-l 

T(m) = - 2A+1 

k=m 

In particular, the mean time to move from one boundary 
state to the other is 


N-l k 

m = E E 

k =0 2=0 


(N — k + e)(fc_j) (■i + l)(fc_j) 
(N — k)(k-i- t-i) (* + e)fe-i+1 


( 6 ) 


The above expression is computationally simpler than 
the product of two hypergeometric functions and involves 
only simple, finite arithmetics. Its expansion in the first 
two powers of e gives (see Mathematical Details): 

T(0) = i + 2 + (7) 

Figure |T| shows the remarkable accuracy of this formula 
for Ne £ [0,1] and N < 100, i.e. the relevant range 
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Figure 1: (Color online) Switching time as a function of e 
for three different values of N. Empty symbols: Numerical 
simulation by a Gillespie algorithm over 10'paths ; filled sym¬ 
bols: numerical solution of the linear system (I5l4ll : Solid lines: 
theoretical expression (0. 


where bi-stability can be observed. The analysis can be 
extended to compute the linear term in e in equation dTJ) 
(see section BA) 

Equations m have been obtained by setting r = 1, 
i.e. by scaling time and e by the factor r. Restoring the 
non-scaled time (f —>■ t/r, e —>■ e/r), we have 

f£ S \ 0) = -T e/r ( 0) 
r 

and in particular, the leading terms of the development 
are 



Figure 2: (Color online) Exact result for the first first passage 

T> 

time (solid line, black) as a function of N for t D = 1/500, 
r = 1 and its comparison to the two solutions provided by 
Biancalani et aid [121]. Figure 5) : dotted curve, blue for e —> 0 

■D 

; dotted curve, red for Ne D > 0.5. 


III. SOLVING THE MASTER EQUATION. 

The mean first passage is one tool to study the stochas¬ 
tic system described by the transition rates m . A com¬ 
plete description can be obtained by solving directly the 
master equation governing the probabilities P(n , t) to ob¬ 
serve n individuals in state A at time t: 



( 0 ) 


1 

e 


2 N — 1 
r N 


-o(-) 

r r 


Therefore, it is possible in principle, by measuring the 
switching time for different system size N, to measure 
independently the parameters e and r. 

Note that the rate coefficients used by Biancalani et 
al. are given in terms of proportions, i.e. r n = N z r 
and e® = Ne. Figure [2] shows the comparison between 
our exact result and the Biancalani et al. approximate 
solutions when this scaling is taken into account, for the 
full range of Ne. It can be observed that the two solutions 
obtained by Biancalani et al. and their overlap can be 
recovered from the exact solution we provide here. 

In a yet unpublished article, Saito and Kaneko 1!)| 
have also computed the switching time for this stochastic 
system. Their method consists in obtaining an approxi¬ 
mation for the residence time toj in each state j begin¬ 
ning from state 0 and then summing up these residence 
times to obtain the switching time. Their analytical re¬ 
sult for the switching time has a very different form that 
the relation ([6]) and doesn’t seem amenable to easy com¬ 
putation of the interesting limiting case JVe C 1. How¬ 
ever, their formula produces the same numerical results 
than the relation © of this article. 


dP{n , t) 
dt 


W + (n — 1 )P{n — 1, t) — W + (n)P(n : t ) 
W~ (n + l)P(n + 1, t) — W~ (n)P(n, tj8) 


We note that the above stochastic system does not need 
a moment closure approximation, i.e. the equation for 
the fcth moment involves only moments of order lower 
than k. Therefore, a hierarchical system of equations can 
be established to derive all the moments of this system. 
The probability generating function is a powerful tool to 
investigate such Master equations fl8l . [20| . The PGF is 
defined as 


N 

<f{z,f) = {z n ) = Y J p Mz n 

n —0 


and contains the most complete information we can have 
on the given stochastic process: all the moments and 
probabilities can be obtained from its derivatives at ei¬ 
ther z = 1 or z = 0. The equation governing the PGF 
can be extracted from the master equation © (see sec- 
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Figure 3: (Color online) The PGF function 0(z, t) as a func¬ 
tion of z at times t € {0, 1, 2,4, 8,16, 32, 64,128, 256}/(128e) 
for N = 100 and e = 0.01. Solid lines: theoretical expression 
cn». Circles: solution obtained by the numerical resolution 
of the Master equation © and computation of its PGF. 



Figure 4: (Color online) The stationary probabilities P(n, oo) 
as a function of n for N = 100 and various e. Solid lines: exact 
expression im, symbols: numerical resolution of the Master 
equation, e = 0.01 (blue circles), 0.1 (green squares), 1 (red 
diamonds), 2 (diamonds, cyan) and 4 (x, purple). 


tion IV Cl) and reads: 


A. Stationary probabilities. 


<90 

~dt 


-z(z- 1) 


2 d 2 <ft 
dz 2 


(z- 1) [(N-l-e)z-(N 
eN(z - 1)0 


1 + 0] 


90 

dz 


( 9 ) 


The solution of equation © can be exactly computed (see 
section IV Cl) as the superposition of polynomial eigen¬ 
functions 


N 

<j>{z,t) = Y C n (j) n (z)e Xnt 

n =0 


( 10 ) 


where the eigenvalues are 

X n = —n(n — 1 + 2e), 
the eigenfunctions are polynomials in z 

N 

0n(O = E ^ - Z t 

k—n 


The stationary probabilities attained at large times are 


p ^={ N r)% { ^r <U) 

(see section IV Cl) and their comparison to numerical solu¬ 
tion of the Master equation is displayed in figure 01 Note 
the qualitative change of behavior at e = 1. Expression 
(fTIT) is eq uivalent to the expression found by Biancalani 
et al. [12| in the continuous approximation, with the ad¬ 
vantage of being well defined for all n, including n = 0, N. 
In particular, for eN <C 1, 


P{n, oo) 


(1 - fr^-ic)/2 + 0(e 2 ) 


Ne 

2n(N—n) 


0(e 2 ) 


n = 0,N 
n ^ 0,7V 


where H m is the harmonic number 


and the coefficients C n depend on the initial condition. 
The initial condition we use here is the same as in the 
previous section, i.e. P(n, 0) = S Uy q which implies that 
0(z, 0) = 1. The exact expression for the coefficients 
flfe , C n and their product are given in the section IV Cl 
The agreement between the solution ([TUI) and the direct 
numerical solution of the Master equation is displayed in 
figure 01 

The PGF contains the most complete information on 
the stochastic process under investigation. Some quanti¬ 
ties of interest extracted from it are given below. 


B. Factorial moments. 

For the purposes of experimental measurements of the 
parameters, other dynamical quantities can be of inter¬ 
est. The most robust of these quantities are the factorial 
moments 

{{n, <?)) = (n(n - 1 )...(n -q + 1)) 

where (n, q) is used to denote the decreasing Pochham- 
mer symbol. The factorial moments are obtained by sue- 
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cessive derivation of the PGF 


((n,q)) =q\ 


dzQ 


Z = 1 


(-i)vE^ 


a q e Xit 


2=0 


( 12 ) 


Note that the gth factorial moment involves only q + 1 
eigenfunctions. The two first factorial moments are 


*> ^ t(‘- 


(n(n- 1)) = 


N(N — 1) 
2 

1 + e 


l + 2e 


— e 


„-2(l+2e)t 


1 + 2e 



For Ne -C 1, only the two first terms in the sum (IT21) 
contribute significantly to the factorial moments for t > 
1. In particular, for large times, 

((n,q)) (iV, c?)-- 


Figure 5: (Color online) T eq as a function of e for different 
values of N. Solid lines: theoretical expression CD ; symbols: 
numerical resolution of the master equation (blue circles N = 
100; green squares N = 50 ; red triangles N = 25). Inset: 
comparison between the exact expression (1 1411 (solid lines) 
and its approximation (1 15 II (dashed lines) for Ne < 1 and N = 
100, 50 and 25. 


C. Equilibrium time. 

Finally, we can define an equilibrium time T eq by 
studying the dynamics of the decrease in P(0,f) or in¬ 
crease in P(N,t). The measure we choose to use here 
is 

pOO 

T eq = {P(N,oo)-P(N,t)}dt (13) 
Jo 

which is a generalization of the mean first passage time 
(see IV Cl ). The expressions for the two boundary proba¬ 
bilities are found to be 

N 

P(0,t) = 

n=0 

N 

P(N,t) = 

n= 0 

and therefore 

N 

T eq = (-fEW A " ( 14 ) 

n= 1 

For Ne < 1, eq. CD is approximated by 

Teq = Te ~ 4 ( i?Ar - 1 _ 2 + n) 

Figure [5] displays T eq as a function of e and its comparison 
to numerical solution of the master equation. 


IV. CONCLUSION. 

As discussed in the introduction, noise induced bi¬ 
stability has been intensely investigated, specially in ge¬ 
netic networks. In general, the chemical Master equa¬ 
tions are too complex to be solved exactly and various 
approximation techniques have been developed to tackle 
this problem. In some cases, exact analytical solutions 
have been obtained using the probability generating func¬ 
tion. Shahrezaei and Swain [2l|| have studied a three 
stage model of simple gene expression (DNA state, RNA, 
Protein) and obtained the protein number distribution. 
Grima et al. [zj | have investigated the steady state dis¬ 
tribution of a two component (DNA state, Protein) ge¬ 
netic feedback loop and have been able to obtain exact 
analytical results using the PGF technique. In the first 
case, the PGF equation is a first order partial differential 
equation and can be solved by the method of character¬ 
istics. In the second case, the model can be reduced to 
two coupled one component systems and the PGF equa¬ 
tion reduced to two ordinary coupled first order differ¬ 
ential equations. Chemical Master equations analogous 
to these cases could in principle be investigated with the 
same technique. 

In this work, we have extended the investigation by 
Biancalani et al. mi of another noise induced bistable 
system which belongs to the second class of models dis¬ 
cussed above. First, we have obtained the exact solution 
for the mean first passage time which is the main result 
of the above cited article. Second, we have solved the 
full master equation associated with this system and ob¬ 
tained other useful quantities for parameter estimations 
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of such systems. We have obtained these results for the 
original, discrete system without resorting to the Tay¬ 
lor expansion of the Master equation in powers of 1/TV. 
Discrete solutions have the advantage of being clearly de¬ 
fined and avoid spurious effect happening at the bound¬ 
aries, specially for the interesting case of small e. More¬ 
over, these solutions involve only simple arithmetic and 
are easily computed. 


V. MATHEMATICAL DETAILS. 

A. Series expansion of the exact solution of the 
switching time. 


B. Solution of Biancalani et al. for the switching 
time. 


In non scaled time, the Biancalani et al. solution is 


T ns (0) = 


X 


1 2N<2 f 1 N 6 ' 3 - 1 ^ 

r'l + 2e'/r' 2 1 \2’ r’ ’ 2’ 1 + 2e'/r’) 


2-Fl 


(I TV—- -■--— 

\2 r' ’ 2’ 1 + 2e'/ 


r' 


where the rates e' and r' are related to the rates e,r used 
in this article through: 

e 1 = Ne ; r' = N 2 r 


The exact solution © contains a double sum, where 
only the terms i = 0 contain e^ 1 factors. Separating 
these two contributions, the solution becomes: 

fYnl = — (!)fc (TV —fc + e)fc 

^ TVe j— (1 + e)k (N~k) k 

+ y' {N — k + ( i + l)(fe_j) 

1^1 ( N ~ k )(k-i+ 1 ) (* + e )k-i+1 

Expanding the first sum to the first order in e necessitates 
only simple expansion in factors of the form to/ (m + e) = 
1 — e/m + 0(e 2 ) and leads to 


C. Deriving and solving the PGF equation. 

PGF. The equation for the evolution of the PGF is 
obtained by multiplying the master equation © by z n 
and summing over n |2.li |. This operation leads to 

^ = ((z n+1 - z n )W+(n)) + <(^- 1 - z n )W~(n)) 

(16) 

The rates W ± (n) are polynomials of second degree in n 
and by the definition of the PGF, 



1 

e 


Hn -i + 2 


N — 1 
TV 


where the Harmonic number H m = Y™ i (1/t). Evaluat¬ 
ing the second sum for e = 0 results in 


N-l k 

i(N — i) 1 

fc= l i= l v ' 


Adding the two contributions results in (eq0: 


T( 0) = - + 2 


TV-1 

TV 


Application of the above rule to equation m leads to 
equation (0. 

Eigenfunctions. Equation m can be transformed 
into a hypergeometric equation by a change of variable 
x = (z — 1) _1 . It is however much simpler to use the 
fact that by definition, the function 4>(z,t ) is a polyno¬ 
mial of degree TV in z and search for the eigenfunctions 
of equation (0 in term of polynomials of the following 
form: 

N 

cfn{z) = Y J <^-z) k 

k -0 


The next term in the series expansion of T(0) is found to 
be 

Note that algorithmically, the computation of T(0) (ex¬ 
pression ©) necessitates only the calculation of TV ratios 
of the form (m+ l)/(m + e) and (m + e )/to which can be 
stored in an array. The T(0) involves then only multipli¬ 
cations and sums of these elements. The Hypergeometric 
function on the other hand is defined as 


OO 

2 F 1 (a,b- c;z) = ^2 

n =0 


( a )(n) (P') (n) Z? 
(c)(n) n\ 


and its efficient implementation requires specific algo¬ 
rithms. 


4>(z,t) = C n (t> n (z)e Xnt 

n—0 

Insertion of these polynomials into equation (0 shows 
that non-trivial solutions (i.e. ^ 0) are possible only for 
the eigenvalues 

X n = —n(n — 1 + 2e) n = 0,1,..., TV 


which leads to a one term recurrence relation on the co¬ 
efficients a% : 


0 (k < n) 
1 


a k +1 = 


(TV — k)(k + e) 


(k + l)(fc + 2e) — n(n — 1 + 2e) 


Ofc (n < k < N) 
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As it can be noticed, (f> n is written as polynomial in pow¬ 
ers of (1 — z) and not z. This choice is not arbitrary: it is 
this change of variable which allows to obtain a one term 
recurrence relation between the coefficients a%. Writing 
c f> n as a polynomial in z leads to a two terms recurrence 
relation which is much more intricate to solve exactly. 

The coefficients can be computed in explicit forms: 


= (-) 


k—n 


N — n\ (e + n)( fe _ n ) 


k-nj (2e + 2n) (fc _ r 


(n < k < N ) 


( 17 ) 

Alternatively, the eigenfunctions can also be given in 
terms of the hypergeometric function: 

(j) n {z) = (1 - z) n 2 -Pi (n - N, n + e; 2n + 2e; 1 - z) (18) 

The amplitudes C n depend on the initial condition. For 
P(n, 0) = S n ,o and therefore (f>(z, 0) = 1, the amplitudes 
obey the triangular linear system 

Co = 1 
k 

= 0 (fc>0) 

n =0 


which can be explicitly solved 


Cn = 


N 

n 


( e )(n) 


and therefore, 


Cnal = (-) 


k—n 


N' 


(2e + n — l)(n) 

(e)(fc) 2e + 2n — 1 


(19) 


nj (2e + n)(fc) 2e + n — 1 


Stationary probabilities. As all eigenvalues except Ao 
are negative, for large times the PGF is simply 

4>{z) = 2 Fi{-N , e; 2e; 1 - z) 

where we have used the hypergeometric representation 
(eq. fl~8l) of the eigenfunctions. Using the relations 


2 Fi(—to, b; c; 1) = 


( c ~ k)(rr 

( c )(m) 


we recover the relation m on the stationary probabili¬ 
ties. 

Factorial moments. Using the above expression, the 
factorial moments are 


((n, 


= w <?) E (-) 1 y ; 

^ \ij (2e + i) (q) 2e + i-l 


(e)(ij) 2e + 2i — 1 


,Xi t 


Equilibrium times. Many different measures can be 
used for the equilibrium time of the system. The expres¬ 
sion we use 


pOO 

T eq = (P(N, oo) — P(N, t)) 
Jo 


dt 


( 21 ) 


is the extension of the mean time to absorption to the 
case when the boundary state is not absorbing. The rea¬ 
son is the following: If the state N were the only absorb¬ 
ing state, whatever the initial condition to, P(N,t) —>• 1 
as t —> oo. The probability of survival until time T, be¬ 
ginning in the state to is 

Q(m, T) = 1- P{N, T) 

and the probability density of not being absorbed during 
[T, T+dt] is therefore — drQ{m, T ). Therefore, the mean 
time to absorption is 


T(to) = — 


poo 

/ Td T Q{m,T)dT 

J o 

poo 

/ (1 -P(N,T))dT 

Jo 

p oo 

/ {P(N, oo) — P(N, T)) dT 

Jo 


dF „ , , n i a )(n){b)(n) „ , , \ 

— 2 F 1 (a,b;c;z) = —— - 2 Pi(a + n, b + n; c + n; z) 

CLZ \C)(n) 

we obtain 


p, , 1 <T4’ 

PM = 5 dF 


We see that in the case of an absorbing state N, our 
definition of T eq and the mean time to absorption are 
the same. We continue to use T eq as a measure of the 
equilibrium time when N is not absorbing. 

Probabilities. The probabilities are extracted from 
the PGF by collecting the coefficients of powers of z: 


N 


P(n, t) = E K exp(A k t) 


k —0 


z =0 


As 


f N)( n ) ( e )(n) ( e )(JV-n) , , 

j n! (2e) (n) (2e + n) ( Ar- ra) 
(2e)(n.)(2 e + n)(jv-n) = (2 e)jv 


where 


N 


^H-mE 


3 ) a k 
n I 0 
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